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Simple models play a key role in the 
microstructural anal^-su; of mechanical 
failure in compoMtcs and other materi- 
als with complex and often disordered 
micrustructures. Although equal luad- 
sharing-modeh arc amenable to rigor- 
ous statistical analysis, problems with 
local load enhancements near failed re- 
gions of the material have so far re- 
sisted exact analysis. Here we show for 
the first time, that some of the simpler 
of these local-load-sharing modek can 
be solved exactly using a sub-stochastic 
matrix method. For diluted fiber bun- 
dles with local load sharing, it is possi- 
ble to find a compact expression for 
the characteristic equation of the sub- 
stochastic matrix, and from it obtain an 
asymptotic expansion for the largest ei- 
genvalue of the matrix. This in turn 
gives the asymptotic behavior of the 



size effect and statistics of the fiber- 
bundle models. Wc summarize these 
results, and show that the important 
features of the exact result can be ob- 
tained from a single scaling analysis we 
had developed previously. Wc also 
compare the statistics of fracture with 
Ihe appropriate limiting extreme-value 
survival distribution (a Gumbel distribu- 
tion), and, as previously indicated by 
Harlow and Phoenix, note that the 
Gumbel distribution performs quite 
poorly in this problem. Wc comment 
on the physical origin of this dis- 
crepancy. 
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1. Introduction 



It has been known since the pioneering work of 
Chaplin [1] and well known since the classic work 
of Griffith [2] that randomly occurring flaws or 
weak links effectively determine the obseived ten- 
sile strength of materials. Early on it was realized 
that the dependence of failure upon the weakest 
part of a material structure gives rise to non-Gaus- 
sian statistics for fracture stress and strain. These 
developments lead to the classical period of the de- 
velopment of the statistics of extremes by mathe- 
maticians such as Dodd [3], Frechet [4], Fisher and 
Tippett [5], von Mises [6], Gnedenko [7], and 
Gumbel [8], 



Following the work of Duxbury et al. [9-1 Ij, 
there have been many attempts to use random net- 
work models to determine the statistics and size 
dependence of material breakdown [12-16]. These 
calculations have in many cases elucidated the gen- 
eral behavior and size dependence of breakdown, 
but few exact results have been produced. 

Perhaps the simplest model that shows the statis- 
tics of brittle failure has been the pure-flaw, chain- 
of-bundles model of Harlow and Phoenbc [17] 
which has been studied by Harlow [18] and more 
recently by Harlow and Phoenix [19]. In this model 
there is a series or chain of m structurally and 



337 



Volume S>9, Number 4, July-August 1994 

Journal of Research of the National Institute of Standards and Technology 



statistically independent bundles of n elements 
each as shown in Fig. 1, where the vertically ap- 
plied uniaxial stress is shared by the surviving verti- 
cal fibers (bonds). Each element or fiber is 
independently present with probability p and ab- 
sent with probability /=!—/?. The survival proba- 
bility of the chain of bundles is then the survival 
probability of a single bundle raised to the power 
m. The main difficulty in this analysis is calculating 
the survival probability of a single bundle. The ex- 
tension of this theory to the survival of two-dimen- 
sional networks is straightforward and amounts to 
the approximation that cracks or flaws only exist 
and break along the direction transverse to the di- 
rection of the applied stress. Following Harlow and 
Phoenix, we assume the local-load-sharing model 
for a flaw (crack) of length n to be 



Ot,p = <T(,(l-i-n/2), 



(1) 



which is to say that the entire force applied to the 
cluster is concentrated at the tip (on the fibers ad- 
jacent to each end of the flaw or vacancy cluster). 
Failure of any surviving bond (all of which have the 
same strength) leads to a rip which causes failure 
of the entire bundle. Solution of this mode! re- 
quires finding the bond (weakest link) which expe- 
riences the largest stress enhancement and that 
stress which would break this most stressed bond. 



Fig. 1. A one dimensional array of intact bonds (fibers) and 
flaws (vacancies). The tensile stress «■, is applied vertically, 

Duxbury, L£ath, and Beale [11] showed how a 
one-dimensional model such as this could be used 
as a simple model for fracture or breakdown of a 
two-dimensional network. If one considers that 
cracks or flaws only exist and break horizontally, 
then the two-dimensional model becomes that il- 
lustrated in Fig. 2 (i.e., no horizontal bonds break). 
Then if we impose spiral boundary conditions, 
where the last site in a row is connected to the first 
site of the previous row on the other side of the 
sample, then the N X-N network problem is re- 
duced to a one-dimensional chain of N^ fibers (or 
bonds) in parallel like that in Fig. 1. 
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Fig. 2. A two dimensional lattice with only horizontal cracks. 
Spiral boundary conditions iUenti^ each site gn the right edge 
with ihc site on the left edge of the previous row. 



2. Single-Cluster Calculation: The Sub- 
Stochastic Transition-Matrix Method 

In this calculation one assumes that the weakest 
link is the pure flaw or cluster of vacancies of the 
largest size that exists in the sample. The survival 
probability is then closely related to the probability 
that in a sample of length L{ = N^} that there are 
no clusters of vacant bonds (flaws) of size greater 
than some prescribed value n. Using a generating- 
function technique Duxbury, L-eath, and Beale [H] 
calculated exactly the asymptotic form of the prob- 
ability to be 



c,X")H^-pr7 



(2) 



in the limit of large L, It is now possible to rederive 
this result while introducing the sub-stochastic 
transition-matrix method. Following Harlow [18] 
we define all possible endings of a fiber bundle of 
length L + 1 and the way in which those endings 
may be generated from a bundle of length L and 
the probabilities of those endings. Since there are 
no allowed clusters larger than n , the allowed bun- 
dle endings or distinctive endings at a particular 
site are an occupied site (1) followed by O^r^n 
vacancies (0) so that these distinct endings 
are spanned by the basis vector <^l^ = 
i4\t)><kioj,4>o«n)-—,4>(w..si))L where the last element 
contains n zeros. 
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Then the probability increment for going from a 
cluster of size r to a cluster of size r' on the next 
site is included in the matrix product 



A/^ = 



p p p ... p 


<^i) 




<^i) 


/ ... 


0(10) 




<^10) 


/ •.. •.. : 
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: i •-. 
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; 


... / 


<^(ll>..4)) 




<^JU...C 



j; u 



(3a) 



; + i 



which is the same as the matrix A/ operating; times 
on the probability vector <^ for the starting site, or 



A/«^-M^^'</po. 



(3b) 



The probability CL{n) that there are only clusters 
up to size n in the entire bundle (or network) of 
size L is thus 



Q(fl)=X(c/.,)L, 



(4) 



where the sum is over all the elements of i^l. The 
simplest and most natural boundary conditions are 
periodic ones where C;.(«) becomes the trace 



C/.(n)=/r(W^)=XA*^ 

k 



(5) 



since the 1st and Lth sites must be the same, where 
Kk are the eigenvalues of A/. We find the eigenval- 
ues of M via its characteristic equation 



D„=detlA//A-/| = 



{a 


-1) 
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... 
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a 
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-1 





... 
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■*^ 


-I 










4*> 
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-1 



= 0, 



(6) 



where a =pJX and 6 =//A =(1 — /j)/A. A cofactor 
expansion of the determinant £>„ about its last row, 
gives immediately the recursion relation 



D, = -£>„„, + (-l)"ai7". 



(7a) 



where D„-i is the n xh determinant for clusters up 
to size (n — 1). With Da = (a — \), the solution, upon 
iteration of Eq. (7a), is the characteristic equation 

(-ly D„=(a -\)+ab +ab^ + -+ab" = 0. (7b) 

Summing this geometric series we obtain 

A"*'-X"*'+pf*'=0. (8) 

This equation is the characteristic equation of M 



times (A -/), so there is an additional spurious 
root at/ (since we are interested in the largest root, 
this does not affect the analysis). Since M is primi- 
tive and non-negative, its largest eigenvalue is real 
and distinct and it can easily be seen that all the 
eigenvalues are less than 1 and Aro„ approaches 1 
for large n . Thus we set Ama = 1 - c and expand Eq. 
(8) to lowest order in e and f, which gives us 



\^^l-pf*' + 0(f'"). 



Then, we obtain 



Q(;i)sA;J;„s(l^pf*y, 



(9) 



(10) 



which confirms the result Eq. (2) by the sub- 
stochastic transition-matix method. 

In order to find the failure probability as a func- 
tion of applied stress, we use the load-sharing rule 
Eq. (1), coupled with the fact that the failure of the 
bond carrying the largest local stress nucleates 
catastrophic failure, and thereby use the relation 



f7\J(r=^l+—, 



(") 



where en, is the breaking strength of a single fiber. 
Note that we could have used a variety of other 
load-sharing rules here, and for example the same 
expression with n raised to an arbitrary power is 
also of physical interest. This result combined with 
Eq. (10) yields the probability S{a-) that a fiber 
bundle will survive at stress a 



S((T) = i\-pf'^-'f 



(12) 



For large n and L, this becomes the modified 
Gumbel form, introduced previously [10,11] in the 
analysis of the random fuse network. Although 
Ci.(n) in Eq. (10) becomes a Gumbel distribution 
in n, the substitution of «(o-) from Eq. (11) pro- 
duces a modified Gumbel distribution that is signif- 
icantly different from a Gumbel form in a in the 
high-reliability tail of the distribution. This modifi- 
cation is di-scussed further in Sec. 4. 



3. Double-Cluster Calculation 

Several authors [12,16,19,20] have suggested that 
the most critical defect is not a single cluster of n 
vacancies, but rather a double cluster (double co- 
linear crack) of n vacancies separated by a single 
occupied site located anywhere within the rt -I- 1 



339 



Volume 99, Number 4, July-August 1994 

Journal of Research of the National Institute of Standards and Technology 



adjacent sites. Such a double crack is shown in Fig. 
3a. This candidate for the most critical crack is ap- 
pealing because the stress enhancement at the in- 
terior occupied site grows as n in network models 
rather than as /i"^ (as for the edges of a single 
crack in a two-dimensional network) and because 
the increased entropy of the (n + 1) locations of the 
occupied site makes it more probable. Thus, fol- 
lowing Harlow and Phoenix [19], we consider the 
probability of bundles of length L with repeated 
double cracks (and single cracks when the occupied 
site is at either end) not exceeding n vacant sites in 
any two adjacent cracks separated by a single site. 
These repeated double cracks are shown in Fig. 3b. 



O) 



b) 



Fig. 3. ■) A double duster of €/i total vacancies plus one i.so- 
lated occupied bond, b) Repeated, civerlapping double dusters; 
each pair of clusters a<> indicated by the brackets contains €n 
total vacancies plus one isolated bond. 

Harlow [18] showed that this problem is 
amenable to analysis by the sub-stochastic transi- 
tion-matrix method. The load-sharing rule is still 
given by Eq. (1) as before but now n is the sum of 
the number of vacant sites immediately on the left 
and right of any isolated intact bond or fiber. Thus 
it is necessary to keep track of not only the number 
of vacancies in the cluster being considered but 
also those in the previous vacant cluster. There are 
now (n + l)(n -l-2)/2 distinct endings that must be 
considered at a site (or bundle ending); these are 
given by the basis vector <^ = <^)),^(io)...<^ia..o); 
<^ 101 ).<^iaio)...<^ioia.,o);<^{iooi)...<^iooio...n);...;c^io.x)))' where 
there are no more than n total vacancies in any 
element. With this ordering of states the « =4 sub- 
stochastic matrix for this problem, for example, is 
given by Eq. (13). 

With periodic boundary conditions Eq. (5) still 
holds and we again analyse the largest eigenvalue 
of Af«. As a technical point, note that since we are 
using periodic boundary conditions, we can always 
start the matrix process at a surviving bond, and so 
the endings considered above include all possible 
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(13) 



survival configurations (we don't have to consider 
configurations which start with O's). 

A great simplification in the characteristic equa- 
tion 



M,<^ = A<^, 



(14) 



where tf> ^^^ the eigenvectors of M„, is possible 
since most of the rows of Mn contain only a single 
non-zero element /. This gives, for example 
f4>i\) = X(kio) or &0()) = (f/A)^i) = ^(i()). By such rela- 
tions, we can eliminate all the rows of M„ except 
the rows with/>'s corresponding to the reduced ba- 
sis vector ^ = <f>(j),<^ioi),<^iooi),.. .,(^100,^)1). The result- 
ing (n+i) equations give an (n+l)x(n + l) 
matrix M', which satisfies the reduced characteris- 
tics equation 
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where a =plX, and b =//A. This M' matrix can be 
considered as a new transition-matrix which adds a 
cluster at a time rather than a bond or fiber at a 
time. Thus we have the characteristic determinant 
equation. 



D„=det 
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For small n these determinants can be evaluated 
directly. For example, 



Do~ -l+a,Di = l-a -a^b and 
D2= - 1 +a +ab +a^b^-a^b\ 



(17) 



But by expanding the determinants Eq. (16) by 
rows and columns, we can show that there is an 
exact recursion relation 

bD„ (fl J> ) =sD„-i{abJb) - D„^i{ab^J) ), (18) 

where Dm{ab'Jb) is D„{ajb) with a replaced by 
{ab% and where 



s^l+b-a^-" 



(19) 



Note that s=s{ajb/i)—s(abjbjn-2)=s{ah'^jb/t 
-4) which is key in the solvability of the recursion 
relation Eq. (18). After some detailed analysis, we 
have found (see [21] for details), that this recur- 
rence equation may be solved. The resulting char- 
acteristic equation is given by 

for n 53 and odd. While for n even S4, we find, 

y n/2-l y nl2-\ 1 



In Eqs. (20a) and (20b), 



2*.-=- 



^± Vj^-4^ 



(20c) 



The key quantity s is given by Eq. (19) above. 
Equations (20a,b) are the exact expressions for the 
characteristic equation of the original M in Eqs. 
(13) and (14). 

Again we find that the largest eigenvalue of M is 
near 1 for large «, so we set A = 1 - c and expand 
Eqs. (20a,b) and find that in both cases, to lowest 
order in e and/". 



e = [{n+T)p^-pT*'+0(f"% 



(21) 



Comparing this double-cluster result to the single- 
cluster result Eq. (9) we find the expected (" +2) 



from the possible locations of the single bond in a 
double cluster of size {n + 1). The ( -p) is a correc- 
tion to properly handle the single-cluster cases as 
well as the double-cluster case, since these are in- 
cluded whenever the isolated bond is located at 
one end of the double cluster. It is only important 
for small n . 

In order to check and better understand the 
asymptotic form Eq. (21) of Abu, and the impor- 
tance of the other eigenvalues A, we have made sev- 
eral numerical evaluations of the various 
equations. First, we have numerically found the 
largest eigenvalue Am« of the original, full, sub- 
stochastic transition-matrix M as given by Eq. (13). 
Using the unit vector as a starting vector we re- 
peatedly apply the matrix M to it. Since the largest 
eigenvalue is unique, this process converges expo- 
nentially to the largest eigenvalue. We found in 
general that convergence occurred to six significant 
figures with at most 50 matrix products (even for 
matrices M of dimension (n + \){n +2)/2 = 10,000. 
The sparsity of M with this iterative procedure 
eliminated matrix-storage problems. The results of 
this iterative procedure for Clin) versus n are 
shown in Fig. 4 as solid lines, and the dots give Amn 
with AmM as given by the asymptotic form Eq. (21). 
Good agreement is seen for allp, with a small devi- 
ation in the p =0.2 data. However, a more strin- 
gent test is needed in the high-reliability (large n) 
tail. Thus, in Fig. 5 we plot the quantity 
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Hg. 4. A plot of the double-cluster probability Ci{n ) vs (n + 1), 
forp -0.2, O.S, and 0.8. The dotted line is the asymptotic form 
given by Eq. (24). The solid lines are found from evaluating the 
largest eigenvalue numerically, and by using Ci(n)'-' Ai„. 
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largest eigenvalue is dominant. A numerical test of 
this convergence is shown in Fig. 6, for L = 1000, 
and shows that for large lattice sizes, the most im- 
portant corrections to Ciin) are the higher order 
contributions to \^^, which are of 0(f^^), rather 
than the neglected smaller eigenvalues of M , which 
are relatively unimportant here. 
Finally, we obtain the asymptotic form 

Cl(")={1 -[(/I +2)p^-pr*'+0(f'"'')f. (24) 

Thus, following the same arguments as at the end 
of Sec. 2, we find that to leading order the survival 
probability of the entire network or chain-of-bun- 
dles is 




Fig. 5. Plots of e^* ' W n + 1, for a) p = 0.5 and b)p =0.2. The 
solid lines are the asymptotic form as given by Eq. (22); the 
circles arc the exact values of (1 - A^)//" * ' as obtained by iter- 
ating M numerically. 

The solid straight line versus (n + 1) is the asymp- 
totic result, which is linear in n and this is com- 
pared with the iterated numerical values (circular 
symbols), for p =0.5 and 0.2. For large n, in all 
cases the two calculations agree. But for small val- 
ues of /7<0.5 there appears a minimum in 
(1 — Amai)^*' versus n which corresponds to higher 
order terms inf. In particular the next order term 
in A^„ is 0(f"'^) which would appear as a Oif"^) 
correction in Fig. 5. 

Finally, we test for the accuracy of neglecting ail 
but the maximum eigenvalue Am^. It is possible to 
directly, numerically evaluate the trace of A/'- by 
iteration since it only requires storage of the matrix 
and one vector at any time. Using this method, we 
have evaluated the quantity 



(l-(rr(A/^))"^)/r*'. 



(23) 



su<r)=[i~i^p'-pr^'^f. 



(25) 



which is again of the form of a modified Gumbel 
distribution with slightly different coefficients from 
Eq. (12). 
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Eq. (22). 



which should converge to {l~\mJ)/f*' when the 
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4. Extreme-Statistical Form 

For large L, we can easily find the limiting form 
of C'i.(n) as given by Eqs. (10) and (24) respec- 
tively, 



Ct (/I ) - exp[ -Lp exp( -n log(l//))], (26) 

/.-.SB 



and 



Ci(n) ~ exp[-L((n +2)p^-p)exp(- 



n !og(l//))]. 
(27) 



as the upper limit behavior for large n which is im- 
portant in the high-surviveability tail of the distri- 
bution. This is a Gumbel distribution, as is 
expected from the exponential (/>/" or «/j/"] be- 
havior for the probability of large clusters [Castillo, 
1988]. 

On the other hand, the sur\'iveabiltties 5(c7-), and 
5 (ct), in the h'mit L-*<», are found to be modified 
Gumbel forms [from Eqs. (12) and (25)], 

Sticr) ~ exp[-L;j exp(-(^-l)log(l//))] (28a) 



and, 



S'da) ~ exp[-L(^p^-/7)exp(-(^-l) 

/,-nJO 



log(l//))]- 



(28b) 



That the dominant behavior of 5i(o-), and S\,{o) is 
exp(-Zw4 exp(-f)) as o- tends to zero is essential 
to ensuring that the survival distributions have the 
proper limiting approach to one when the applied 
stress approaches zero. Harlow and Phoenix [19] 
have numerically shown that this high-reliability 
tail can not be well described by a Gumbel form for 
u (such as exp(-Z^ exp(B((7-oti))))- But this 
failure is obvious since the Gumbel form doesn't 
approach one until cr^-— =o, so at sufficiently small 
stresses it must be inaccurate. Nevertheless the 
standard texts on extreme distributions (see, for ex- 
ample [8] or [22]) seem to suggest that the Gumbel 
distribution is the appropriate one in such cases. 
The difference clearly is in the form of the nor- 
mally assumed scaling limw.^.[S(o-)]''=5'(flA'cr + fcjv), 
which fits the shift and slope of the limiting func- 
tion 5((r) at its median, but fails near the high-reli- 
ability limit o- = 0. Ft would seem that instead a 
more general scaling form limv^.[5(n(cr))}'" = 
5(o/vn((r) + bw) must be allowed to at.so include the 
proper high-surviveability limit near zero stress. In 
many practical material-failure problems, this mod- 



ification of the Gumbel form is essential in order to 
correctly represent the important high-reliability 
tail. Note that this is not true of the Weibult distri- 
bution, which as well as being a stable limiting ex- 
treme-value distribution, does have the physically 
correct behavior as stress approaches zero. This is 
one good reason why the Weibull distribution is a 
very robust form in the analysis of failure problems. 
We suggest that the family of modified Gumbel 
distributions of the sort Eq. (28), should be simi- 
larly robust, in contrast to the conventional Gum- 
bel distribution which is of limited use in the 
analysis of the statistics of material failure. 



5. General Scaling Behavior 

The size dependence and general form of the 
limiting distribution can usually be found from a 
back-of-the-envelope scaling calculation which we 
introduced previously [9,11]. First, for the single- 
cluster calculation, consider the probability /'£(«) 
of finding a cluster of size n in a sample of size L. 
The order of magnitude of this probability is 



Vtin) ^Lp^rip^Lpf, 



(32) 



since, for normalization, 1p^f'=p, and since 
there are L different locations in the sample where 
a cluster could be located. For the maximum clus- 
ter size to be expected in a sample of size L we set 



Pi.(n)~Lpf'^ = ], 
and obtain the size dependence 



^rnax — 



logL 
Jog(7/) 



(33) 



(34) 



or, from the load-sharing rule Eq. (1), the break- 
down stress oi 



02 



»■}. 



1 



1 



l+ftwJ2 J InL 



(35) 



21n(l//) 



scales to zero logarithmically in the thermodynamic 
limit. A similar argument for double clusters gives 



Fi(n)^Lin + l)p'r, 



(36) 



since there are (« + 1) places to put an isolated 
bond in the n -double-clusters. We then obtain the 
same limiting form Eq. (35) for crc(L), although 
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there are additive {log(logL)) corrections in the 
double-cluster case. The logarithmic scaling law in 
turn implies that the failure statistics is of the dou- 
ble-exponential form given in Eq. (28). The 
Weibull and Frechet distributions always give 
power-law size scaling. These qualitative argu- 
ments are very powerful and are confirmed by the 
rather elaborate, exact calculation described here. 
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